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Abstract 


This paper presents a forviulatwn for the identification of linear 
multivariable systems from single or multiple sets of input -output data.. 
The system input-output relationship is expressed in terms of an ob- 
server , which is made asymptotically stable, by an embedded eigenvalue 
assignment procedure. The prescribed eigenvalues for the observer may 
be real, complex, mixed real and complex, or zero. In this formulation , 
the Markov parameters of the observer are identified from input- output 
data. The Markov parameters of the actual system are then recovered 
from those of the observer and, used to obtain a. stale space, model of the 
system by standard realization techniques. The basic mathematical for- 
mulation is derived, and extensive numerical examples using simulated 
noise- free data are presented to illustrate the proposed method. 


Nomenclature 


Symbol 

Dimensions 

Description of symbol 

,4 

n x n 

system matrix (discrete) 

A 

n x ii 

observer system matrix (discrete) 

A„ 

n x 11 

system matrix in observable canonical form 

Aur 1 ) 


polynomial of delay operators associated with 
output 

n 

ii x in 

system input influence matrix (discrete) 

B 

ii x (in -f- q) 

observer input influence matrix 

B' 

n x in 

partition of D associated with input u ( i ) 

Ki 

n x m 

partition of D associated with input u(i) for a 
( lea< 1 beat, ( ) 1 >ser ver 

B(q l ) 


polynomial of delay operators associated with 
input 

B* 

n x in 

transformed partition B f of B for a MI MO 
system, B* - TB f 

b 

n x 1 

system input influence vector for a SISO system 

b 

n x 2 

observer input influence vector for a SISO systei: 

bo 

n x 1 

system input influence vector for a SISO system 
in observable canonical form 

y 

71 x 1 

partition of b associated with scalar input 

b* 

n x 1 

transformed partition // of b for a SISO system, 
b * - Tbf 

<>; 

scalar 

ith element of the transformed vector b* 

vi; 

1 x rn 

ith row of the matrix B* 



c 

q x n 

system output matrix 

c* 

q x n 

transformed system output matrix for a MIMO 
system, C* = CT l 

c 

1 X 71 

system output vector for a SISO system 

Co 

1 x n 

system output vector for a SISO system in 
observable canonical form 

c* 

1 x n 

transformed output vector for SISO systems, 
c* = cT~ ] 

<■* 

scalar 

ith element of the transformed output vector c* 

r* 

C (i) 

q x 1 

ith column of the transformed system output 
matrix C* 

D 

q x m 

direct transmission matrix 

d 

scalar 

direct transmission term for SISO systems 

Em 

srn x m 

a matrix of identities and null matrices, 



x m Omx(s- 1 ) rn 

E q 

rq x q 

a matrix of identities and null matrices, 



ts 

II 

1 

X 

0 

X 

1 

i 

H(r) 

qr x ms 

an r x s block data matrix of Markov parameters 
for realization 

I 


identity matrix 

M 

n x q 

observer gain 

M d 

n x q 

deadbeat observer gain for a MIMO system 

M* 

n x q 

transformed observer gain for a MIMO system, 
M* = TM 

injT 

(») 

1 x q 

ith row of the matrix M* 

m 

scalar or 

number of inputs or 


n x 1 

observer gain for a SISO system 

m a 

n x 1 

observer gain for a SISO system in observable 
canonical form 

mi 

n x 1 

deadbeat observer gain for a SISO system in 
observable canonical form 

rnf 

scalar 

ith element of observer gain in observable 

canonical form 

m* 

n x 1 

transformed observer gain for a SISO system, 
77i* —Trn 

m i 

scalar 

ith element of the transformed vector m* 

n 

scalar 

order or assumed order of a linear system 
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U c 

scalar 

number of prescribed complex conjugate pairs of 
eigenvalues 

n r 

scalar 

number of prescribed real eigenvalues 

0 


null matrix 

p 

scalar 

number of observer Markov parameters, C 13. 
CAD , . . C A B , also referred to as window 

width in the identification algorithm 

Pi 

scalar 

coefficients of characteristic equation for a SISO 
system 

<1 

scalar 

number of outputs 

(l~ l 


one-time step delay operator 

/?", R m 


space of 77-, r/-, and 777 -dimensional real- valued 
vectors 

r 

scalar 

pole radius of prescribed eigenvalues in the 
complex plane 

T 

n x n 

a similarity transformation matrix for ^4. .4 — 

r 1 at 

U.V 

rq x n, ,s*m x n 

orthonormal matrices obtained from the singular 
value decomposition of the Hankel matrix //(C)) 

u(i) 

m x 1 

input to system at time step i 

u(i - p) 

rrip x 1 

p-time step input history vector 

u(i — n ) 

run x 1 

77 - time step input history vector 

v{i) 

(m + q) x 1 

vector containing input and output served as 
“input” to observer at time step / 

,:(*) 

7 7 , x 1 

system state vector at time step i 

x(i) 

7 7, x 1 

estimated system state vec tor at time step i 

x(i) 

77 X 1 

state estimation error vector at time step i 

Y T 

q x m 

system Markov parameter, Y- = CA r D 

Y r 

q x (777 + f/) 

observer Markov parameter, Y T = C A T D 

v(U 

1 T 

q x 777 

partition of Y r associated with input u(i) 

yW 
1 r 

r/ x q 

partition of Y r associated with output y{i) 

V(i) 

q x 1 

output of system at time step 7 

y(i ) 

g x 1 

estimated output at time step i 

y(i - n) 

77 ? x 1 

n-tiine step output history vector 

y{i - p) 

qp x 1 

p-timc step output history vector 

a 

77 x 1 

SISO observer parameter vector associated with 
input for real eigenvalue assignment 


3 



OL C 


<*d 


Om 


(Y 


(Y r 


<Ad 


ft 


ftr 


0a 


ftll 


ft 

ftr 

(]d 

ftrn 

r r (i-i) 

rd(i-i) 
r m (i - 1) 
r(i- i) 
rr(i-i) 


n x 1 
n x 1 
n x 1 
q x /mi 
<7 x 7/m 
(/ x ///.//, 

q X 7/77/ 

77, X 1 
7/ X 1 

n x 1 
71 x 1 
q X c/7/ 
x qn 
q X CJ71 
q x gn 
(2n+ 1) x 1 
(2n + 1) x 1 
(2 n+ 1) x 1 
(2 n + 1) x 1 

(rn + 9)7/ +- 771 x 1 
(777 + q)n + rn x 1 


counterpart of a for SISO complex eigenvalue 
assignment 

counterpart of a for SISO deadbeat eigenvalue 
assignment 

counterpart of a for SISO mixed eigenvalue 
assignment 

MIMO observer parameter matrix associated 
with input for real eigenvalue assignment 

counterpart of a for MIMO complex eigenvalue 
assignment 

counterpart of a for MIMO deadbeat eigenvalue 
assignment 

counterpart of a for MIMO mixed eigenvalue 
assignment 

SISO observer parameter vector associated with 
output for real eigenvalue assignment 

counterpart of ft for SISO complex eigenvalue 
assignment 

counterpart of ft for SISO deadbeat eigenvalue 
assignment 

counterpart of ft for SISO mixed eigenvalue 
assignment 

MIMO observer parameter matrix associated 
with output for real eigenvalue assignment 

counterpart of ft for MIMO complex eigenvalue 
assignment 

counterpart of ft for MIMO deadbeat eigenvalue 
assignment 

counterpart of ft for MIMO mixed eigenvalue 
assignment 

combined vector of (j)(i — 1), <p(i — 1), and u(i) for 
SISO real eigenvalue assignment 

counterpart of r(i — 1) for SISO complex eigen- 
value assignment 

counterpart of T(i — 1) for SISO deadbeat 
eigenvalue assignment 

counterpart of F(?' — 1) for SISO mixed eigenvalue 
assignment 

combined vector of <j>{i — 1), <p{i — 1), and u(i) for 
MIMO real eigenvalue assignment 

counterpart of F(i — 1) for MIMO complex 
eigenvalue assignment 
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r d (i - 1) 


r m (i — i) 
7 
7 r 


I'd 


fm 


7(i) 


>(*) 


7tf(0 

7m(0 

I 
7 c 


7,/ 


7 m 

7(0 

7c (0 

7 d(i) 

7m (0 

A 

Ac 

Am 


(m + r/jn + m x 1 
(rn + q)n + rn x 1 
(27) + 1) x 1 
(277 -(- 1) x 1 
(277 + 1) X 1 
(277 + 1) X 1 
(277 + 1) X 1 

(2n + 1 ) x 1 
(277 + 1) x 1 
(2?) + 1) x 1 
q x ( 77 ? + q)n + m 
q x (777 + <7)77 + rn 
q x (777 + 77)77 + m 
r/ x (m + < 7 ) 7 ? + rn 
q x (m + qr)n + m 

q x (m -f </)n + rn 

q x (m + q)n + rn 

q x (m + r/)r?, + rn 

n x n 
n x n 

n x n 


counterpart of T(i - 1) for MIMO deadbeat 
eigenvalue assignment 

counterpart of £ (i - 1 ) for MIMO mixed eigen- 
value assignment 

vector of combined observer parameters for SISO 
real eigenvalue assignment 

counterpart of 7 for SISO complex eigenvalue 
assignment 

counterpart of 7 for SISO deadbeat eigenvalue 
assignment 

counterpart of 7 for SISO mixed eigenvalue 
assignment 

vector of combined observer parameters esti- 
mated at time step i for SISO real eigenvalue 
assignment 

counterpart of 7 (/) for SISO complex eigenvalue 
assignment 

counterpart of 7 (i) for SISO deadbeat eigenvalue 
assignment 

counterpart- of 7(7) for SISO mixed eigenvalue 
assignment 

matrix of combined observer parameters for 
MIMO real eigenvalue assignment 

counterpart of 7 for MIMO complex eigenvalue 
assignment 

counterpart of 7 for MIMO deadbeat eigenvalue 
assignment 

counterpart of 7 for MIMO mixed eigenvalue 
assignment 

matrix of combined observer parameters esti- 
mated at time step i for MIMO real eigenvalue 
assignment 

counterpart of 7 ( i ) for MIMO complex eigenvalue 
assignment 

counterpart of 7(7) for MIMO deadbeat eigen- 
value assignment 

counterpart of 7 (i) for MIMO mixed eigenvalue 
assignment 

diagonal matrix of prescribed real eigenvalues 

block diagonal matrix of prescribed complex 
eigenvalues 

counterpart of A for mixed real and complex 
eigenvalues 
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prescribed real or complex eigenvalues 


A, 


X( T ) 

n x 1 

A, (r) 

n x 1 

a!? 

n x 1 

A (r) 

-t.m 

m x rn 

A (r) 

qxq 

x (T) 

Am 

run x m 

\ {r) 

A<\rn 

run x rn 

A (r) 

Am. in 

run x rn 

A (r) 

Aq 

qn x q 

x ( t) 

A(\q 

qn x q 

\ {t) 

Arn.q 

qn x q 

E 

n x n 


scalar 

d T) , ^ r) 

scalars 

* {T) 
— ijri 

rn x rn 


qxq 

<p( l - 1) 

n x 1 

(pc(i - 1) 

n x 1 


vector of powers of prescribed real eigenvalues for 
a SISO system 

counterpart of A^ for SISO prescribed complex 
eigenvalues 

counterpart of A^ for SISO prescribed mixed 
eigenvalues 

diagonal matrix of r- powers of real eigenvalue A ? ; 
repeated m times 

diagonal matrix of r - powers of real eigenvalue A, 
repeated q times 

matrix of powers of prescribed real eigenvalues 
associated with input for a MIMO system 

counterpart of A^/ associated with input for 
MIMO complex eigenvalue assignment 

(t) 

counterpart of A ?n associated with input for 
MIMO mixed eigenvalue assignment 

matrix of powers of prescribed real eigenvalues 
associated with output for a MIMO system 

counterpart of A^ associated with output for 
MIMO complex eigenvalue assignment 

counterpart of A q associated with output for 
MIMO mixed eigenvalue assignment 

diagonal matrix of positive singular values 

real part of a complex eigenvalue 


elements of the matrix 


" », w 

i 



1 

3 

i 

1 

' ’ H 

* -H 

l 


"* U/ j 



( 7 ” ) 

diagonal matrix of a- repeated m times 

diagonal matrix of a- repeated q times 

transformed p-time step input history vector 
<p(i — 1) = 3u(z — p) for SISO real eigenvalue 
assignment 

counterpart of <j>{i — 1) for SISO complex eigen- 
value assignment 
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< fim(i ~ 1) 

n x 1 

counterpart of <j)(i — 1) for SISO mixed eigenvalue 
assignment 

<j)(i - 1) 

nrn x 1 

transformed p-time step input history vector 
(j){i — 1) — % n u(i - p) for MIMO real eigenvalue 
assignment 

<M* - l) 

nrn x 1 

counterpart of <f>(i — 1) for MIMO complex 
eigenvalue assignment 

~ 1) 

nrn x 1 

counterpart of 0(z — 1) for MIMO mixed eigen- 
value assignment 

<p(i - !) 

n x 1 

transformed p-time step output history vector 
< p(i — 1) — Qp(i — p) for SISO real eigenvalue 
assignment 

<Pc.(i ~ 1) 

n x 1 

counterpart of i p(i — 1) for SISO complex eigen- 
value assignment 

'Pm (t 1) 

n x 1 

counterpart of p(i — 1) for SISO mixed eigenvalue 
assignment 

p(i ~ ! ) 

ng x 1 

transformed p-time step output history vector 
ip(i — 1) = — p) for MIMO real eigenvalue 

assignment 

fS ~ !) 

nq x 1 

counterpart of (p(i — 1) for MIMO complex 
eigenvalue assignment 

1 pm{i 1) 

nq x 1 

counterpart of p(i — 1) for MIMO mixed eigen- 
value assignment 

LU'i 

scalar 

complex part of a complex eigenvalue 

J T) 

—i,m 

m x rn 

diagonal matrix of u;- 7 ^ repeated in times 

J T) 

qxq 

diagonal matrix of ud repeated q times 

3 

n x p 

Vandermonde- like matrix of real eigenvalues for 
SISO systems 

% 

n x p 

counterpart of 3 for SISO complex eigenvalue 
assignment 

%n 

n x p 

counterpart of S for SISO mixed eigenvalue 
assignment 

—m 

ran x mp 

Vandermonde-like matrix of real eigenvalues 
associated with input for a MIMO system 

— cjn 

mn x mp 

counterpart of for MIMO complex eigenvalue 

assignment 

— 771,771 

mn x mp 

counterpart of for MIMO mixed eigenvalue 

assignment 

c> 

-q 

qn x qp 

Vandermonde-like matrix of real eigenvalues 
associated with output for a MIMO system 

> 

—c,q 

qn x qp 

counterpart of ( 3 q for MIMO complex eigenvalue 
assignment 
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%nji V n x < IP 

counterpart of c 3 q for MIMO mixed eigenvalue 
assignment 

■(?(?') (m + q)n + rn 

x (m + q)n + rn 

projection or variance matrix for recursive least- 
squares estimation 

Abbreviations: 


ARMA 

auto-regressive moving average 

ERA 

Eigensystem Realization Algorithm 

MIMO 

multiple-input multiple-output 

SISO 

single-input single-output 


Introduction 


The aim of learning identification is to provide methods to improve identification of the 
system model as additional information about the system becomes available. The techniques 
are in the time domain, and the system information comes in the form of input-output data 
from either multiple experiments or a single experiment of extended duration. Originally, 
the idea of learning identification was motivated by the fact that for system identification of 
flexible structures, multiple experiments are often performed with the hope that the averaged 
data can reduce the effects of irregularities such as measurement noises, repetitive disturbances, 
and slight nonlinearities. This motivates the development of learning identification to improve 
identification results effectively from multiple experiments. An early technique for identification 
of parameters from multiple experiments was formulated in reference 1. Learning identification 
in the present form identifies the Markov parameters from general input-output data (ref. 2). The 
Markov parameters are then used to obtain a state space model of the system by a realization 
procedure, e.g., the Eigensystem Realization Algorithm (ERA) (refs. 3 and 4). The learning 
identification procedures presented in reference 2 require input-output data from a large number 
of experiments of generally short duration. The procedures identify as many Markov parameters 
as the number of data samples in each experiment, and the number of data samples that can be 
used is constrained by the number of Markov parameters desired to be identified. In practice, 
there may be substantially more data samples in each experiment than the number of desired 
Markov parameters. Therefore, these techniques are not efficient in the sense that they do not 
necessarily make use of all available input-output data. This motivates the development of 
identification algorithms from a single set of input-output data of extended duration. Learning 
identification is closely related in concept and technique to learning control, where the motivation 
is to develop control laws that improve tracking error based on repeated execution of a task 
(refs. 5 11). 

An identification procedure from a single set of input-output data is developed in reference 12 
by means of an auto-regressive moving average (ARMA) description of the original system 
in state space format via an observer. An important distinguishing feature of the approach 
presented in reference 12 as opposed to previous development is that the system is identified 
indirectly by an observer, which is made asymptotically stable by an eigenvalue assignment 
procedure. The discrete- time eigenvalues are required to be real, distinct, with magnitudes less 
than one. The recursive formulation of reference 12 extends the repetition domain concept used 
in learning control and identification to shifting time intervals. It is based on procedures that 
identify system Markov parameters for indirect learning control and repetitive control (refs. 9 
and 10). The use of Markov parameters in system identification is discussed in reference 13. In 
this paper, the identification technique is generalized to allow assignment of complex eigenvalues. 
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This generalization is particularly important when the order of the system is large, since it 
permits a more even distribution of asymptotically stable eigenvalues inside the unit circle by 
using the entire complex plane. 

The basic contributions of this paper are as follows: First, a simplified reformulation of the 
original identification algorithm with placement of real eigenvalues is presented. Second, the 
formulation is extended to the case of complex eigenvalue assignment, which is also applicable 
to the general case of assigning both real and complex eigenvalues. Third, a version of this 
identification procedure using a deadbeat observer with poles placed at the origin is formulated. 
This case is of particular interest since it makes use of a minimum number of data samples, and 
the number of identified observer Markov parameters are reduced to a minimum set. Fourth, 
an extensive numerical study is provided to illustrate the basic characteristics of the algorithm. 
The deterministic technique developed here is applicable for data from either a single set or 
multiple sets of experiments. Because of the complexities in the formulation of the identification 
algorithms, the case of single- input single-output systems will be first described. The results 
are then extended to the case of multiple-input multiple-output systems. This paper gives a 
more detailed presentation of the results in reference 14, with additional examples. In order 
to study the exact nature of the identification procedure under ideal circumstances, t his paper 
is confined to purely deterministic results. In the presence of process and measurement noises, 
the relationship between the identification algorithm with a deadbeat observer presented in this 
paper and the stochastic Kalman filter algorithm of reference 15 is established in reference lfi. A 
procedure to improve observer and Kalman filter identification results by whitening the residual 
sequence is presented in reference 17. Often of interest in practice is the identification of a 
model in a prescribed frequency range. Such a development of the algorithm is formulated in 
reference 18. 

The general outline of the paper is as follows. The procedure with real eigenvalue assignment, 
which is first presented in reference 12, is reformulated here using a modified mathematical 
formulation. The modified formulation allows direct extension of the procedure to the case 
with complex eigenvalue assignment. A special case of the identification procedure when all 
eigenvalues are placed at the origin is then presented. For clarity, the formulation for single- 
input single-output systems is presented in the main body of the paper, except for the case of 
mixed real and complex eigenvalues assignment, which is presented in appendix A. Extensions 
of the identification procedure to multiple-input multiple-output systems arc parallel to the 
developments for the single-input single-output case. The multivariable case is presented 
in appendix B. The truss structure used in the numerical example section is described in 
appendix C. 


Mathematical Preliminaries 

The following general mathematical formulation is applicable to both single-input single- 
output (SISO) and multiple-input, multiple-output (MIMO) systems. This section introduces 
the basic concepts and establishes some general mathematical relations, which are used to derive 
the identification algorithm in subsequent sections. 


System Description 

In this section, the relationship between the state space model and a particular auto- regressive 
moving average (ARMA) model of a linear system is presented. This relationship is particularly 
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useful for developing an identification procedure. First, consider a general discrete multivariable 
linear system expressed in the state space format 

x(i + 1) — Ax(i) + Bu(i) 
y(i ) = Cx(i) + Du(i) 

where x(z) G R n \y(i) G R q ,u(i) G R m . Let x(0) denote the initial state at i = 0. An input- 
output description of the above system can be obtained from equation (1) as 

y(i) = CA‘x( 0) + ^2 CA i ~ T ~ ] Bu(t) + Du{i) (2) 

r=0 

Note that the first term on the right-hand side of the above equation is dependent on the initial 
condition ;r(0). The products CA l ~ r ~ l B denoted by Yi- r - 1 , together with D , are known as 
the Markov parameters of the system. From equation (2), the input-output description of the 
system with zero initial conditions becomes 

i - 1 

y(i) = ^ Y r u(i - r — 1) + Du{i) (3) 

T—0 

where y(i) is expressed in terms of Yq up to Y^\ and the direct transmission term D. In general, 
this description requires i + 1 Markov parameters to describe the output at time step i. If the 
system is asymptotically stable such that the Markov parameters Y p , Y p + i, Y p + 2? ••• can be 
neglected for some p , then at time steps i > p, the input-output description can be approximated 
with a finite set of Markov parameters as 

p - 1 

y(i) « ^ Y r u(i — r — 1) + Du(i) (4) 

r= 0 

It is important to note that for a finite dimensional system, there is only a finite number of 
independent system Markov parameters. Therefore, the system Markov parameters used in the 
description of equation (4) are not necessarily independent. For sufficiently damped systems, 
equation (4) is a valid description of the input-output relationship provided that p is chosen 
sufficiently large such that the approximation holds. However, for lightly damped systems, such 
as large flexible space structures, the ARMA model would require a very large number of Markov 
parameters, which would not be computationally attractive for system identification. In fact, 
if the system is unstable or marginally stable, such a description is no longer possible. In the 
following, a procedure is described to express the state space model in equation (1) as an ARMA 
model with a finite number of Markov parameters. The Markov parameters can be shown to be 
those of an observer system that is made asymptotically stable by eigenvalue assignment. This 
observer model is then used to develop an identification method for the system described by 
equation (1). 

To construct an observer model, add and subtract the term My(i) to the right-hand side of 
the state equation in equation (1) to yield 

x(i + 1) = Ax(i) + Bu(i) + My(i) — My(i) 

= (A + MC)x(i) + (B + MD)u(i) - My(i) (5) 
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Define 


A = A + MC 
B = [B + AID , -M] 

> 

u(i) 

v(i) = 

2/(0 J 

J 

Then the original system becomes 


( 6 ) 


x(i + 1) = Ax(i) + Bv(i) 
y(i) = Cx(i) + Du{i) 


The input-output description of the system with zero initial conditions is 


( 7 ) 


i-l 

V(i) = 5Z Y i—T— l , -( T ) + Du(i) 


r=() 


( 8 ) 


where 

y /: - r _ 1 =CA i - r ~ l B 

If the system is made asymptotically^stable by placing of the eigenvalues of the matrix ~4 
such that the Markov parameters Y p , F ?H _ l5 F />+2 , can be neglected for some p, then at 
time steps i > p, the input-output description can be approximated with a reduced set of p + 1 
Markov parameters [Y q, V j, . .., Y p ^\, D}. The following equality 


p - 1 

*/(0 = X! Y r h '{i - r - 1) + Du(i) (7 > p) (9) 

7"=0 

then approximately holds. If the original system is observable, then for any system matrix ,4, it 
is always possible to find a matrix M such that the desired eigenvalues of A are placed in any 
particular (symmetric) configuration. For the case of lightly damped systems, this procedure 
can transform the set of an otherwise large number of Markov parameters to an approximately 
equivalent reduced set {Tq, Y . ... Y p _ }, D} by selecting appropriate eigenvalues for ,4. 
Furthermore, for a sufficiently large /?, the influence of a nonzero initial condition on the output 
at time steps i > p can be neglected. The model of equation (9) is used to develop the 
identification method presented herein, and the eigenvalue assignment step is achieved implicitly 
through processing of the measured input-output data. To see that equation (9) is a special auto- 
regressive moving average model, it can be rewritten as 


' r 

2/(0 + C ( A + MC) T My(i -t-1) = ^2 C(A + MC) T (B + MD)u(i - 


r=() 


1)T Du(i) (10) 


Defining a delay operator q 1 applied to a variable z{i) to be q~ x z{i) = z(i - 1), the above 
equation can be written in the usual deterministic ARMA model format 

A (? _1 ) 2/(0 = B ((T 1 ) u{i) (11) 
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with the polynomials of the delay operators A (q *) and B (q given as 
A ( 9 _1 ) = I A CMq~ l A CAM q~ 2 + • • • + CA p ~ l Mq- p 
B (V 1 ) —Da CB'q ~ X A CAB'q ~ 2 A • • • + CA p ~ x B'q- p 

where A — A + MC and B f = B + M D. 

Relations of the System to an Observer Model 

The role of the matrix M in the above development can be interpreted in terms of an observer 

model. Consider the system given in equation (1). It has an observer of the form 

x(i + 1) = Ax(i) + Bu(i) - M[y(i) - yii)} 

y(i) = Cx(i) + Du(i) (12) 

It can be shown from equations (12) and (1) that 

x(i + 1) = Ax(i) + Bu(i) — MC\x(i) — x(i)\ 

= (^4 + MC)x(i) + Bu(i) — M[y(i) — Du(i)\ 

= {A + MC)x(i) + {B + MD)u(i) - My(i) (13) 

Defining the state estimation error x(i) = x(i) - x(i ), the equation that governs x(z) is 

x(i + 1) = Ax(i) + Bu(i) - [(A + MC)x(i) + (B + MD)u{i) - My(i)] 

= (A + MC)x(i) (14) 

If system (1) is observable, then M may be chosen to place the eigenvalues of A + MC in any 
desired (symmetric) configuration. In particular, they will be placed inside the unit circle in the 
complex plane. From equation (14), if M is chosen such that A + MC is asymptotically stable, 

then lim x(i) = 0; i.e., the estimated state x(i) converges to the true state x(i) as i approaches 

i — ►oc 

infinity. Equation (13) then becomes 

x(i + 1) = (A + MC)x(i) + (B + MD)u(i) - My{i) (15) 

which is exactly the same as equation (5). 

From this analysis, matrix M can be interpreted as an observer gain. The parameters 
Y x _ T _l = C~A l ~ T ~ x B in equation (8) are then the Markov parameters of an observer system; 
hence they are now referred to as observer Markov parameters. In the identification process, 
these are the parameters to be identified. Once they are identified, the actual system Markov 
parameters can be recovered. There is an algebraic relationship between the Markov parameters 
of the observer system and those of the actual system. This result is established in the following 
section. 

Relations Between the Markov Parameters of the Observer and the Actual System 

As before, let the Markov parameters of the observer system be denoted by Y r , and the 
Markov parameters of the actual system by Y r . Recall that 
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Yr = CA T B 


= [C(4 + MC) T [B + MD), —C(A + MC) T \1] 
= [Y { r l) ,Y { T 2) } 


(16) 


From the second equation in equation (16), the Markov parameter CB of the system is simply 


Y 0 = CB = C(B + MD) - {CM) I) 


— y^ 1 ) -L v 


y n 
1 0 u 


(17) 


To obtain the Markov parameter CAB , first consider the product 


Hence, 


Y = C(A + MC)(B + MD) 

= CAB + CMCB + C(A + MC)MD 


Yi = CAB 

= Y ^ + rj ) 2 ) Y 0 + y[ 2) d 

Similarly, to obtain the Markov parameter CA 2 B, consider the product y ^ 

Y ( 2 1} = C(A + MC) 2 (B + MD) 

= C{A 2 + MCA + AMC + MCMC)(B + MD) 

= CA 2 B + CMCAB + C{A + MC)MCB + C(A + MC) 2 MD 


Therefore, 

y 2 = CA 2 B 

= Y ( 2 1) - CMCAB - C{A + MC)MCB - C(A + MC) 2 MD 
= y[ l) + Y l, 2) Y 1 + Y < 2) Y 0 + Y (2] D (19) 

By induction, the general relationship between the actual system Markov parameters and the 
observer Markov parameters can be shown to be 


Y t = 




+ g F m 

1=0 


y r _ f _i +y {2) d 


(20) 


For a noise-free finite-dimensional system, knowledge of a sufficient number of actual system 
Markov parameters is adequate to deduce a state space realization of the system of interest. 
Physical aspects of the model such as natural frequencies, damping ratios, and mode shapes can 
then be found. 
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Identification Theory for Single-Input Single-Output Systems 

In the following, an identification method is developed to identify the coefficients of an ARM A 
model that is made asymptotically stable by an embedded eigenvalue assignment procedure. The 
coefficients of the ARMA model are precisely the observer Markov parameters formulated in the 
above section. For simplicity, consider the case of a single-input single-output system in the 
state space format: 


x(i + 1) = Ax{i) + bu(i) 

y(i ) = cx(i) + du(i) (21) 

where ( i ) £ H n . and u{t) and y(i) are scalars. The system matrix A is an n x n matrix, b an 
n x 1 column vector, calx n row vector, and the direct transmission term d is a scalar. The 
input-output description of this system is given as in equation (9): 

P~ 1 

y(i) = ^2 YMi ~ r - 1) + du{i) (22) 

T = 0 


where 

Y T — cA T b — c(A 4- mc) T [b + md , —m] = [cA T b f , —cA r ?n] (23) 

The observer gain m in this case is an n x 1 column vector. Recall that v(i) contains both the 
input u(i) and the output y(i). For i > p, equation (22) can be rewritten as an approximate 
ARMA model: 

p — 1 p-l 

y(i) = y^(cA T b') u{i-T - 1) + du(i) - ^(cA T rn) y (z - r - 1) (24) 

T — 0 T = 0 


Derived in the following section is an algorithm that computes the coefficients c.A^b' and cA T rn 
of the ARMA model in equation (24) and simultaneously places the eigenvalues of A in prescribed 
locations so as to make the ARMA model asymptotically stable. These eigenvalues may be real, 
complex conjugate pairs, a combination of both, or zero (deadbeat). 

SISO Real Eigenvalue Assignment 

This is the simplest case, where all the prescribed eigenvalues are real and distinct. The 
eigenvalue assignment procedure can be derived by noting that for desired real and distinct 
eigenvalues of A, one has for some nonsingular matrix T 

A = T _1 AT (25) 


where A is a diagonal matrix of n prescribed eigenvalues, 


A = 


^2 


(26) 


An 
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For simplicity, the blank spaces denote zero elements. Then the products cA T b f and 
become 

cA T h' = cT ~ 1 A r Tb f cA T m = cT~ 1 A r Tm 

If the elements of c* = cT~^ , b* = Tb f , and m* = Tm are written explicitly as 


c = 


1 ( 2 


b * = 


V 


i 

* * 

i 

*2 

777* — 

77? 2 



i 

2 * 

1 


then the product cA T b in equation (24) may be expressed as 


cA T l/ = c*A T b* = 


.. 


c* fr* 
*■ ir'n 


A5 


\ T 


Similarly, 


cA 


rn 


= -c* A 7 


777 ~ 


T n \ 


With the following simplified notations 



-c n m n 





T 

a = 


8 t = 


A (r > = 


r*h* r*h* 

c l°l C 2°2 


l"‘l ~ i: 2" l 2 


A[ A5 x; t 


equation (24) becomes 


P - 1 P-! 

y(i) = c\ T ^ A ( r \i(i — r — 1) + ^ A ^y(i — r — 1) + du(i) 

r = 0 r— 0 


or simply 


y{i) = 7 r r (i - 1) 


cA T m 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 
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which is in linear form with the unknown parameters in a T , 0 T . and d with 


T 

7 = 


a T 0 T d 


F(i - 1) = 


4 >(i - 1) 
<f(i ~ 1) 
«(*) 


( 34 ) 


where 


P-1 


</>(i — 1) = A^u(i - r — 1) = ^su(i — p) 


T — 0 


P-1 


( 35 ) 


<p(i ~ 1) = A(r) y(* - T - 1) = - p) 

T — 0 

The matrix Q is a Vandermonde-like matrix of prescribed real eigenvalues of magnitudes less 
than 1: 


Aj _1 \ p ~ 2 

1 \P~ 2 

a 2 a 2 


Ar 1 Ar 2 ... a« i 

and the p x 1 input and output history vectors u(i — p) and y(i — p) are defined as 


Ai 1 

A2 1 


( 36 ) 



1 

'a' 

1 

1 


y(i - p) 

u(i — p) = 

u(i — 2) 

y(i -p) = 

y(i - 2) 


1 

1 

3 

1 


y(i - 1) 


( 37 ) 


Note that equation ( 33 ) is in linear form; thus the unknown observer parameter vector 7 
can be solved for directly from input-output data. For on-line computation, however, recursive 
solution is often preferred. Let 7(1) denote the estimated parameter vector at time step i. The 
standard recursive least-squares solution to equation ( 33 ) is simply 


7 (0 = 7(* - 1) + 


R(i-2)r(»- 1) 
i + r(i-i) T 3 ?(f- 2 )r(i-i) 


y(i) - 7 r (i - i)r(» - 1) 


»(» -l) = R(i -2)- 


*(< - 2)r (i - i)r (* - i) t 3 ?(* - 2) 
1 + r(i- l) T 5 R(f - 2)r(i - 1) 




( 38 ) 


with an arbitrary initial guess 7(0), and 9 ft( — 1 ) is any positive definite matrix. Other recursive 
parameter estimation algorithms may be used to replace the standard least squares at this step, 
e.g., the projection or instrumental variable methods (refs. 19 and 20 ). The above algorithm 
identifies the parameter vector 7, which consists of the products c?b*, -c*ra* (i = 1, 2, . . . , n), 
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and d. These products together with the assigned eigenvalues A ? (i = 1, 2, . . . , n) can he used 
to reconstruct the observer Markov parameters Y T (r = 0, 1, 2, . . .) as 

Y r — cA r b = cA T [b f , — m] 



= a T \^ d T X^ (39) 

Finally, the actual system Markov parameters can then be recovered from the above' recon- 
structed observer Markov parameters according to equation (20): 

r r = F< 1) + £Y< 2 V T _,_ 1 + F< 2) D 

7 = 0 

= « T \W + I3 T ^ A (T 'F r _ + \Wd\ (40) 

where = ct 1 A^ T ) and Y ^ 

SISO Complex Eigenvalue Assignment 

With the general mathematical framework developed for real eigenvalue assignment, the 
procedure for complex eigenvalue assignment can be similarly derived by replacing equations (25) 
and (26) with their counterparts for complex conjugate pairs of eigenvalues A; = cq ± jujj 
(i = l, 2, n/2). Namely, A = T~ 1 A C T, and 
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Furthermore, associated with the complex conjugate pair A; = < 7 j ± juj , . write 



LOi 

T 

1 

A 

( r ) 

u i 

-LOi 

J 


(t) 

L“ w i 

T' 

i 


then the product cA T b' in equation (24) becomes 


( 42 ) 


cA T b' = c*A T c b* 


= n\ T] (cp; + dt?) + ^ (cK - + 4 T ‘ lr :V'.'l + <$»!) + 4 T) - cJ*S) 

+ • • • + + 4K) + «-iK - <W,-i) 


= Q^ r) 


(43) 


Similarly, 


/iT * * r * 

cA rn = c A,>m 




(cj m| + {%m%) + u;[ r) {c\m* 2 - c%m\) + a) 2 1 (r^rn^ + cjmj) + (c^m\ - c.pm^) 


(t) 


,( T ) 


I (^") / * * I * , (t) / * 

+ %/ 2 V n-l m n—l + c n m n) + </ 2 { c n 


(t) 

n /2 


-1 m n - r n m n l) 


= -^a 5 t) 


(44) 


where 


cJ6J+c.^ rJ6* 


' c 2 6 1 






a t = - 


A< r) = 


c l”‘I + <- 2 w 2 r l m 2 - c 2 m ! • • • r n \ m n I + c > l n r n-l m n - < tn n-l 

(^’> = 1, uf =o) 


T (r) , ,(r) (t) (r) 


(f ) (r) 

<T n/2 ^ n /2 


(45) 


The elements c *, b*, and m * (i = 1, 2, n) are defined exactly the same way as in 

equations (28), and n is now necessarily even, since all the prescribed eigenvalues must appear 
as complex conjugate pairs. Equation (24) now becomes 


y(i) = a ( T y \[ T ^u{i - t - 1) + flj A [j\{i - r - 1) + du(i) 

T = 0 T=0 

= lAc(i ~ 1 ) 


(46) 
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which is in linear form, with 


4>c{i - 1) 

7, 7 = (*T (i r c (i - 1 ) = <Pc{i- 1 ) 


where 


Mi - 1) = X < T)u (i -T- 1) = - p) 


■(* _ 1 ) = ^2 X <^y(i - r - l) = l 3 r y(i - p) 


The matrix Cs r is a Vandermonde-like matrix of prescribed complex eigenvalues of magnitudes 
less than unity 


, ,(/>-!) , ip ~ 2 ) 


rri 1 


Jp-2) 


a 2 1 


; 2 0 


^ (? ;: 2) ... 

n/2 n/2 

,.(P-1) ,,(/>- 2) 

. w/2 u 'n/2 


n n Pi 1 


w n/2 


and tlie px 1 input and output history vectors u(i—p) and y(i-p ) arc defined as in equations ( 37 ). 
Let 'yc(i) denote the estimated parameter vector at time step i. The recursive least-squares 
solution for the complex eigenvalue case is obtained by simply replacing 7 (i) by >(<). F(i - 1 ) 
by r f; (i — 1 ) in equations (38) with an arbitrary initial guess 7 ,( 0 ) given, and 'ft(-l) is any 
positive definite matrix 5ft(>. Any other recursive algorithm may be used to replace the standard 
least squares at this step. The algorithm identifies the parameter vector 7 ,., which consists of the 
product sums and differences c*_ l b*_ { + c *b* , c*_ j b* - c^b *_ , , c*_ , m*_ , + c*m * , c *_ , m* - c* rn *_ , 
(* ~ 1- ..., n), and d. These identified parameters, together with the assigned conjugate 

pairs of complex eigenvalues^, (* = 1 , 2 , ..., n/ 2 ), can be used to reconstruct the observer 
system Markov parameters Y r (r = 0 , 1 , 2 . . . .) 

Y T = cA T b = cA T [6',-m] 


= a]'xl T) = Y ( r ] Y [ r ] 


Finally, the actual system Markov parameters can then be recovered from the above recon- 
structed observer Markov parameters according to equation (20) in the same way as the real 
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eigenvalue case: 



( 51 ) 


SISO Deadbeat Eigenvalue Assignment 

If all the eigenvalues of the deterministic observer system are placed at the origin, then the 
Markov parameters of the observer system will become identically zero after a finite number of 
time steps. This is a deadbeat observer. Specifically, 

Y t = 0 (r = n, n+1, n + 2, . . .) (52) 


where n is the order of the system. Let rn ( f denote the deadbeat observer gain, the expression 
relating the input-output of the system and the corresponding observer Markov parameter is 
given by 

n — 1 7i — 1 

y(i) = ^2 ( cA T b ') u(i - t - 1) - ( cA r m d ) y{i - r - 1) + du{i) (53) 

7=0 7—0 


The structure of A can be better seen by considering the system given in equation ( 21 ) in 
observable canonical form: 


0 


-a\ 


by 

1 

0 

-CIO 


b2 


1 

-a$ 

b 0 = 

^3 



0 : 





1 — (in 


bn 


(54) 


T 

Let the observer gain in observable canonical form be denoted by m 0 = [mi m 2 m 3 ■■ ■ m n ] . 
The observer system matrix A = A„ + m„c 0 is simply 


0 

1 


0 

1 0 


1 


— a\ + m\ 
-02 + m 2 
—03 + m3 


1 —a n + rn H 


(55) 


For a prescribed set of eigenvalues \ for A, the observer gain m 0 is unique and its elements are 
given by m z = «, - />, , where p,; are the coefficients of the characteristic equation 


(A — Ai) (A — A 2 ) • • • (A — A n ) = A u + P/i A" 1 + • • • + P 2 A + Pi — 0 (56) 


Let m,g = [mf ■■■ mf t } T denote the deadbeat observer gain for the system in observable 

canonical form. In the deadbeat case, the characteristic equation is simply A n -- 0. Hence, p, = 0. 
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m'i = <ij . The observer system matrix A then becomes 

0 

1 0 

A= 1« (57) 

1 

1 0 

In this case, it is convenient to work with the system in observable canonical form directly. The 
Markov parameters of the deadbeat observer system can be computed as 


Yq — c 0 b 0 — c 0 


Y l = (- oA b 0 = 


b 0 + rnid -mi 


b n + rn' ri d —mi 


bn-l+mf^d -mi , 


Y n— ] = c 0 A n 1 b () = 


b\ + rn { [d 


Y n — V 7J+1 — Y n + 2 * 1 ' — 0 


Equation (24) becomes 

77 — 1 


77— 1 


y(i) — X] (bn-T + u(i - r - 1) - ^ - r - 1) + du(i) 

T = i) T~0 

Defining the parameter vectors 

b n + m' 7 [d b n -\ + m d n _ x d ■■■ b\ + m'(d 


T 

a d = 


( 1 ) 


0 T ,= 


y(l) y 
r 0 r 1 


— — m 


yO) 

r n- 1 


n 


y Y ^ 

r 0 Y 1 


-mV 


:( 2 ) 

77-1 


equation (59) can then be written as 

y(i) - ajiu(i - n) + djy(i - n) + du{i) 

= 7 j’r d (*-i) 


(58) 


(59) 


(60) 


(61) 
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where 




d 


r d (i - 1 ) = 


u(i — n) 

y(* - «) 

«(*) 


(62) 


and the n x 1 input and output history vectors u(i - n) and y(i - n) are defined as 


1 

E£ 

1 

i 


1 

"S' 

1 

" (Si 

Pi 

i 


y(i - n) = 


u(i — 2 ) 


y{i - 2 ) 

i 

1 


r 

<SJ . 

1 

t— 1 

1 


(63) 


The recursive solution to equation (61) is obtained by simply replacing 7 (i) by 7 d(i), and T(i- 1) 
by _ 1 ) i u equations (38). The actual system Markov parameters can be recovered according 
to equation ( 20 ) as 


t — 1 


Y r = Y + J2 Y r )Y r-i~ 1 + Y r >d 

i = 0 


(2), 


71— 1 


= b n -T + mS_ T d - 


«-r d ) 




1=0 


n— 1 


— bn—T ^ r^r— i— 1 


(64) 


where T = mf t _ T = 0 , for r = n, n + 1 , 

A particular feature of the deadbeat algorit hm is that the observer system Markov parameters 
are identically zero after a finite number of time steps. The input-output ARMA relation given 
in equation (22) or equation (24) used in deriving the algorithm therefore holds exactly. This is 
different from the previous cases, where by placing real and complex eigenvalues of magnitudes 
less than unity but greater than zero, the ARMA relation only holds approximately. The degree 
to which the approximation holds depends on the choices of prescribed eigenvalues and the 
window width p, i.e., the number of observer Markov parameters retained to maintain a valid 
approximation. In the deadbeat case, however, the approximation becomes exact, the window 
width p is the order of the system, and the identified parameters contain an exact description 
of the system of interest. 


Realization by the Eigensystem Realization Algorithm 

A state space model of the system from the recovered Markov parameters can be obtained 
by the Eigensystem Realization Algorithm (ERA). The algorithm begins with an r x s block 
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data matrix called the Hankel matrix and denoted by H(r) 


Y r 

*r+l • 

■■ y t+s _ i 

Yt+ 1 

^r+2 •' 

T t+ . s . 

^r+r— 1 

Yr+r " 

^r+r-f.s— 


(65) 


The order of the system is determined by the singular value decomposition of H( 0), 

H{ 0) = UT,V T (66) 

where the columns of U and V are orthonormal, E is an nxn diagonal matrix of positive singular 
values, and n is the order of the system. Defining a <\ x rq matrix El. , and an m x srn matrix 

rp i 

Ej U made up of identity and null matrices of the form 


E T - 
- 


Iqxq Oqx(r-l)q 


El = 


1 m x rn O. 


m w mx(s-l)ra 


(67) 


a discrete-time minimal order realization of the system can be shown to be 


A r = Z~ 1/2 U t H( 1)VY,~ 1 / 2 
B r = E V 2 V T E m 


C r = E%u E 1 / 2 


(68) 


This is the basic ERA formulation. To use ERA in the present identification procedure, the 
entries that make up the data matrix given in equation (65) are precisely the recovered system 
Markov parameters Y T (r = 0, 1, 2, . . .). For further details on the algorithm, the readers are 
referred to various references in the literature, e.g., references 3 and 4. 


Computation Steps 

This section reviews the basic steps involved to implement the identification procedure 
developed in this paper. The related equations are identified in each step of the process. 


Step 1 

Assume an order n for the system to be identified. Choose an order p for the ARMA model, 
and select the prescribed eigenvalues of the observer. For the eigenvalue assignment procedures, 
p is normally several times larger than the assumed order of the system, n. Specifically, the value 
of p chosen must be consistent with the prescribed eigenvalues for the observer, as described in 
the following: 

(a) For real eigenvalues, select n real eigenvalues A (i = 1, 2, . . . , n) such that X 1 ’ ss 0. 
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(b) For complex eigenvalues, the eigenvalues must appear in complex conjugate pairs, 
A, = Cj ± jui (i — 1, 2, . . . , n/2) such that. 


<Jj UJ ij 
(Tj 


0 


(c) For a combination of real and complex eigenvalues, the same rules apply. 

(d) For deadbeat observers, however, all eigenvalues are set to be zero, and p is the same as n. 

The identification equations for the deadbeat case have taken this into account. Therefore, 

no explicit specification of the eigenvalues for this case is necessary. 

Note that for asymptotic stability, all prescribed real or complex eigenvalues must have 
magnitudes less than unity. 

Step 2 

Compute the observer parameters. The appropriate recursive equations used for each case 
are outlined as follows. For observers with assigned real eigenvalues, equations (38) are used for 
the SISO case, and equations (B13) are used for the MIMO case. For observers with complex 
eigenvalues, the recursive equations are obtained simply by replacing 7 (i) by %(i), and F(i - 1) 
by T c (i - 1) in equations (38) for the SISO case, and by replacing 7 by 7c, and F(i - 1) by 
r r (i - 1) in equations (B13) for the MIMO case. For observers with mixed real and complex 
eigenvalues, replace 7(?‘) by 7 m (t), and T(i — 1) bv T m (i — 1) in equations (38) for the SISO 
case, and 7 (i) by 7 m (i), and T(i - 1) by r m (i - 1) in equations (B13) for the MIMO case. For 
deadbeat observers, the appropriate recursive equations are obtained by replacing 7 {i) by 7 y(i), 
and T (i - 1) by F ( i(i - 1) in equations (38) for the SISO case, and 7 (i) by ), and T^(i - 1) 
by T d {i — 1) in equations (B13) for the MIMO case. 

Step 3 

Reconstruct the observer Markov parameters from the identified observer parameters. For 
observers with real eigenvalues, equation (39) is used for the SISO case, and equation (B14) 
is used for the MIMO case. For observers with complex eigenvalues, equation (50) and 
equation (B23) are used, respectively. Similarly, for observers with both real and complex 
eigenvalues, equation (All) and equation (B31) are used. For deadbeat observers, however, the 
identified parameters are precisely the observer Markov parameters, and no reconstruction of 
the observer Markov parameters is needed for this case. 

Step 4 

Recover the system Markov parameters from the observer Markov parameters. The general 
equation is given in equation (20), which is then specialized to various cases. For observers 
with real eigenvalues, equation (40) is used for the SISO case, and equation (B15) is used for 
the MIMO case. For observers with complex eigenvalues, equation (51) is used for the SISO 
case, and equation (B24) is used for the MIMO case. For observers with both real and complex 
eigenvalues, equation (A12) and equation (B32) are used, respectively. For deadbeat observers, 
equation (20) directly applies. 

Step 5 

Realize a state space model for the identified system from the recovered system Markov 
parameters in step 4 above. The basic equations for ERA are summarized in equations (65) 
to (68). 
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Numerical Examples 

The theoretical development sections discussed the use of observers and eigenvalue placement 
to recover the system Markov parameters. The Markov parameters are the pulse response 
samples of a linear system. The fundamental idea in the developed identification procedure 
is to identify parameters of an observer rather than those of the actual system. From the 
observer parameters the true system parameters can be recovered. The observer eigenvalues 
or poles determine the observer pulse response decay rate. Formulations where the prescribed 
eigenvalues are real, complex, mix real and complex, and zero (deadbeat) have been presented. 
By making the pulse response of the observer system decay sufficiently fast through the 1 placement 
of its poles, one can truncate the response after a finite number of time steps. Because of the 
different eigenvalue placement procedures, this approximation will result in different convergence 
properties for each respective algorithm. 

To study the numerical properties of the identification procedure, an analytical model of a 
truss structure is used. The lightly damped structure, known as the Mini- Mast (ref. 21) at 
NASA Langley Research Center, is modeled by its first five modes, with frequencies of 0.80, 
0.80, 4.36, 6.10, and 6.16 Hz. A more detailed description of the system under consideration is 
given in appendix C. The outputs correspond to displacement sensors, and the inputs to torque 
actuators. The input-output data are simulated using random inputs for 6 sec. The system is 
discretized at a sampling rate of 33.3 Hz, and an input-output history of 200 points is recorded 
for system identification, which is performed on a Macintosh Ilci computer. The analytical 
model contains five modes, but practically only three of them are controllable and observable 
from any given input-output pair. 

Single-Input Single-Output Examples 

First, for clarity the case of single-input single-output identification is studied. Basic 
characteristics of the identification algorithm can be seem in the SISO case. For this purpose, 
the first input second output pair is used for identification, which results in a system with 
essentially three identifiable structural modes, i.c., a sixth-order system. Results for the cases of 
real, complex, and deadbeat eigenvalue assignments arc presented. The case of mixed real and 
complex eigenvalue assignment is omitted here since its numerical properties may be deduced 
from those of real and complex eigenvalue assignments. 

Consider the case of real eigenvalue assignment. The identification results for this case are 
reported in figures 1(a) 1(d). Figure 1(a) shows the nominal case where six observer poles are 
placed at ±0.2, ±0.3, and ±0.4. Along with the prescribed pole locations is an estimate of the 
number of samples or window width p that it takes for the observer pulse response to decay to a 
negligible value. In this example, the window width is selected to be 40 points wide, i.e., p = 40, 
so that (±0.2)^, (±0.3) p , (±0.4)^ are negligible. The identification procedure starts with an 
initial estimate of the system order, which for the nominal case the assumed order is six, n = 6. 
Even though the model used is of lOth-order, from any input-output pair the effective order of 
the system is only 6. 

The top left plot of figure 1(a) shows convergence histories of the observer parameter values 
calculated from the recursive least-squares solution given in equations (38). The constant values 
correspond to converged parameters. Since the initial parameters are assumed to be zero, to 
start the algorithm, the projection matrix ft(-l) is set to a large value to reflect the degree 
of uncertainty of the initial guess. The plot on the top right shows the square root of the 
diagonal elements of the variance or projection matrix 3?(i) after 160 iterations of the recursive 
least-squares algorithm. In cases where the exact least-squares solution is obtained and no order 
over-specification occurs, the variance matrix approaches zero. In general, the variance matrix 
provides a measure of the freedom in the uniqueness of the identified parameters. Note that when 
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the identified parameters are not all independent because of order over-specification, the large 
variance values do not imply inaccuracies in the parameter estimates. This merely means that for 
the specified order, the identified set of observer parameters is not unique. The recursive least- 
squares solution is driven by the prediction error shown in the second row of figure 1(a). At any 
time step, the prediction error is defined to be the difference between the true output value and 
the predicted output value computed based on the estimated model available at that time step. 
The initial prediction error is large but quickly goes to zero as the observer parameters converge 
to constant values. For the case of real eigenvalue assignment, from the identified observer 
parameters, the observer Markov parameters are recovered by equation (39). The actual system 
Markov parameters are then computed by equation (40). Using the computed pulse response, 
realization of a state-space representation of the system is performed using equations (65)-(68). 
At this step, the initial assumption made about the system order (n — 6) is verified by counting 
the number of nonzero singular values. Shown in the second row of figure 1(a) is a bar chart of 
the normalized singular values which shows six nonzero singular values. 

The top four plots of figure 1(a) are indicators as to how well the parameters are identified. 
The bottom four plots show 7 results comparing the identified state space model and the true 
system model. Included in this group are comparisons of realized and actual pulse responses; 
actual displacement history used in the identification, and its reconstruction using the identified 
model; and the frequency response functions. There are two curves in each plot; the solid curve 
corresponds to actual data and the dashed curve to reconstruction. When an exact model of the 
system is identified, the two sets of curves overlap. 

To study the effect of order under-specification, figure 1(b) shows the results w'hen the 
observer poles are placed on the real axis, as in figure 1(a), but the assumed system order 
is set to n = 2. This is a case where not enough freedom is allowed in the identification 
procedure. The parameter values, shown on the top left of figure 1(b), do not tend to constant 
values as in figure 1(a). Although the variance is small, the prediction error show r s discrepancies 
between the predicted and actual outputs. Realization using the identified parameters results 
in a system of order two, as shown by the singular value plot. When comparing the impulse 
responses, it is clear that the results are in error. So are the reconstructed displacement and 
frequency response functions. In this case, the algorithm attempts to identify a sixth-order 
system by a second-order model. Figure 1(c) shows the results when the assumed system order 
is increased to four, n = 4. Convergence of the parameters is observed, and the corresponding 
variance is small. The prediction error fluctuates about zero. The realized system order is four, 
as depicted in the singular value plot. Comparing the pulse responses shows very small errors. 
However, the frequency response functions show' that the identified system (depicted by the 
dashed curve) missed the mode with the smallest contribution to the system response. This is 
why the reconstructed displacement, when compared with the actual displacement as showm in 
the lower left plot, shows no visible differences. This example suggests a potential application 
of the algorithm for identification of reduced order models. 

To examine the case of order over- specification, figure 1(d) shows the results when the observer 
poles are also placed on the real axis as before, but the order of the system is over specified to 
be n = 10. Results are similar to those shown in figure 1(a), with two important distinctions. 
First, the parameter variances are now substantially larger than those in the previous cases; in 
fact they are an order of magnitude larger than the identified parameters. Second, the realized 
system order is correctly identified to be 6 even though the initial assumed order is 10. Large 
variances are expected when the identified parameters are not all linearly independent. When 
order over-specification occurs, there are more parameters than necessary to identify the system 
exactly. It is important to observe that at the realization step, however, the system and its order 
are identified correctly. For the case of SISO identification, if the assumed order is less than 
or equal to the true order of the system, as shown in figures 1(a) 1(c), the algorithm returns 
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an identified model with the same order as assumed. However, if the assumed order is more 
than the true order, the algorithm returns a model with the correct minimal order, as shown 
in figure 1(d). The identification procedure, as mentioned earlier, places the observer poles at 
prescribed locations. To verify the proper eigenvalue placement, the observer pulse responses 
are used to realize the observer model, and the recovered eigenvalues are found to be identical 
to the prescribed values in all cases. 

The next group of figures (figs. 2(a) 2(d)) presents results when complex poles are prescribed. 
Six complex poles are placed on a circle with a radius r = 0.5 in the complex plane corresponding 
to the same damping level. The window width is selected to be p = 40. The top left 
plot shows the parameter convergence histories. The identified parameters are now given by 
equation (45) instead of equations (31). The overall performance given in terms of prediction 
error, reconstructed response, pulse response, and frequency response functions is similar to that 
of the real case. Results for the complex case with order under-specification, n = 2, are shown 
in figure 2(b). When the assumed order is increased to four, figure 2(c) shows that the identified 
solution misses the weakest mode of the system. Again, this is consistent with previous results. 
In the complex case when the system order is over specified, n = 10, some of the parameters 
do not converge, as shown in figure 2(d). Nevertheless, the system and its order are identified 
correctly. This points out that there are linearly dependent parameters that are being identified. 
This is also indicated by the large variances computed. 

To study the effect of truncation error when the pulse responses have not decayed to zero 
in the allowed window width, the pole radius in the complex case is increased to 0.9 while 
maintaining the same window width p = 40. Results in figure 3(a) show the parameter values 
drifting, while the variance is relatively small. The correct system order is used iit this example. 
The prediction error is large, and the realization procedure identifies a fourth-order system. The' 
identified pulse response, the reconstructed output, and the frequency response functions are 
significantly different from those of the actual system. The situation can be easily corrected by 
increasing the window width to p = 80 to reduce the truncation error. This is verified by the 
results presented in figure 3(b). 

To eliminate the truncation error, the observer poles can all be placed at the origin. This 
is known as the deadbeat case because the observer pulse responses will go to zero in exactly 
a finite number of time steps. No estimate of the window width is needed, because once an 
assumption about the system order is made, the window width is automatically fixed. Results 
for the deadbeat case assuming the correct order are shown in figure 4(a). These' results are 
similar to the real and complex cases, although the identified parameters are different. In all the 
cases discussed, the same input-output time histories are used for identification. Figures 4(b) 
and 4(c) show the deadbeat case when the assumed order is two and four, respectively. 
Figure 4(d) depicts results for order over-specification. 

As with any numerical method, proper conditioning of the data is important. When 
identifying systems where the magnitudes of the input and output values arc orders of magnitude 
apart, because of the use of different units for example, proper scaling of the numerical values 
is critical. This is true even for simple systems. The results shown here are scaled such that the 
input arid output values have comparable magnitudes prior to application of the algorit hm. 

Multiple-Input Multiple-Output Examples 

Identification of MIMO systems proceeds similarly to the SISO case. The model is the same, 
but now two inputs and two outputs are used for identification. The system is now of order 10. 
As in the SISO case, 200 data points are used in the identification algorithm. 

First, consider the case where all prescribed observer poles are real. Initially the assumed 
order is set to n — 4 with corresponding pole locations at ±0.2 and ±0.3, and the window 
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width p is set to 40. The top row of figure 5(a) shows the parameter convergence histories and 
the variance distribution for the least-squares solution. The parameters seem to have reached 
constant values, but some variations arc still observed. Since the true order is 10, this is a 
case where the assumed order is less than the true order. This lack of freedom in the identified 
parameters prevents the prediction error from converging to zero, as shown in the second row of 
figure 5(a). Counting the number of nonzero singular values, the identified system order is found 
to be eight. The system pulse response, the reconstructed output, and the frequency response 
functions are in error. Figure 5(b) shows the results when the assumed order is increased to 
six. The parameters converge to constant values but the variances are large. Large variances 
indicate that, some redundancy in the identified parameters has occurred. However, this does not 
affect the final answers. The prediction error converges to zero, and the system order is correctly 
identified to be 10. This is in contrast with the SISO case, where the identified system order does 
not exceed the assumed system order. It is important to note that for a given set of poles, the 
pole placement problem generally contains an infinite number of solutions for a multiple-output 
system. This results in additional freedom in the algorithm that is not present in the SISO case. 
It is this freedom that allows the identification of a system with a higher dimension than initially 
sought. One interesting aspect of the MIMO case is that when the observer pulse responses arc 
realized to verify the prescribed pole locations, the apparent observer order is equal to the 
assumed order times the number of outputs. The resulting observer poles are those prescribed 
initially, but they arc repeated as many times as the number of system outputs. A comparison 
of the pulse response, the reconstructed response, and the frequency response functions for the 
second output shows excellent agreement. Results for the first output are similar and not shown 
here. 

Figure 6(a) shows results when the prescribed observer poles are complex. The poles are 
distributed evenly in the complex plane on a circle with radius r = 0.5. The assumed order 
is four. As in the real case, the window width p is set to 40. Results are similar to those of 
the real case in figure 5(a). Figure 6(b) shows the complex case when the assumed order is set 
to six. For the deadbeat algorithm, figures 7(a) and 7(b) show the identification results with 
assumed orders of four and six, respectively, when all prescribed poles are placed at the origin. 
Performance of the identification algorithm is similar to the previously discussed examples. As 
in the SISO case, if the assumed order is higher than the true order, the system can still be 
correctly identified, and the algorithm returns an identified model of minimal order. 

Concluding Remarks 

This paper formulates an algorithm for identification of linear multivariable systems from 
general input-output data. Data from either single or multiple sets of experiments can be used 
to identify or update the system model. For each data set, the initial condition may be arbitrary 
and need not be known. The procedure identifies the Markov parameters of an observer system 
instead of those of the actual system. The actual system Markov parameters are recovered from 
the observer Markov parameters and then used to realize a minimal state space model of the 
system. The embedded eigenvalue assignment procedure is used to specify the observer with 
asymptotically stable poles. The prescribed poles may be real, complex, or mixed real and 
complex. When all the prescribed poles are placed at the origin, this results in an identification 
algorithm with a deadbeat observer. Expressed in linear form, the observer Markov parameters 
can be solved for in one step for off-line computation, or recursively for on-line computation. 
The standard least-squares algorithm, which is used in one step of this identification procedure, 
may be replaced by other recursive parameter estimation algorithms. Identification procedures 
for both single-input single-output and multiple-input multiple-output systems are formulated, 
and numerical examples using noise-free simulated data are presented to illustrate the basic 
characteristics of the developed method. 
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Appendix A 

SISO Mixed Real and Complex Eigenvalue Assignment 

The identification procedure with mixed real and complex eigenvalue assignment can he 
derived in the same way as the case with real or complex eigenvalues. Among n prescribed 
eigenvalues, let n r denote the number of prescribed real eigenvalues A, (i — 1 , 2, .... n r ) and 
n r thenumber of prescribed complex eigenvalues, A, = a r ± juij (i = 1 , 2, ..., n r / 2). Then 
write A = T _ * A m T, and 
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The elements c*, 6*, and m* (i = 1, 2, . . . , n) are defined exactly the same way as ii 
equations (28). Equation (24) now becomes 
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which is again in linear form with 
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The matrix T m is a Vandermonde-like matrix of n prescribed real arid complex eigenvalues of 
magnitudes less than unity: 
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and the p x 1 input and output history vectors u(i — p) and y(i — p) are defined in equations (38). 
The standard recursive least-squares solution for the mixed real and complex eigenvalue case 
is obtained by simply replacing the estimated parameter vector 7 (i) by 7 m (i), and T (i — 1 ) by 
T m (i - 1) in equations (38). The observer Markov parameters Y r (r — 0, 1, 2 . • • ■) can be 
reconstructed according to 
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Finally, the actual system Markov parameters can then be recovered as 
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Appendix B 

Generalization to Multiple-Input Multiple-Output Systems 

The developed identification theory for single-input single-output systems can be extended 
to the multivariable case. Consider the multivariable system in equations (1). The input-output 
relation in terms of the Markov parameters of an observer system is given in equation (9), which 
can be rewritten in ARM A model form as 
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A recursive algorithm that computes the matrix coefficients of the ARMA model, and at the 
same time places the eigenvalues of A at prescribed locations, is derived in the following sections. 
These eigenvalues again may be real, complex, a combination of both, or zero (deadbeat). 


MIMO Real Eigenvalue Assignment 

Let the prescribed eigenvalues of A = T~^AT be denoted by (? = 1, 2, • • * , n). Then the 
products CA r B f and CA r M become 

CA r B f = CT~ l A r TB f CA T M - CT~ l A r TM (B2) 

If the elements of C * = CT” 1 , = TB f , arid M * = TM are written explicitly as 






b ([) 


m *T 

m (D 



r 

-1 


b m 


*T 

m (2) 


c* = 

r * r * 
C (l) (2) 

' * C/ V 

(n) 

B* = 

M* = 

(B3) 

f> (n) 

rn* ( r 

(«) 


where denotes the itli column vector of the matrix C*, and and m {i = 1, 2, • • • , n) 
denote the fth row vectors of the matrices B* and M * , respectively, then the products in 
equation (Bl) may be expressed as 
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where 4,m and A^ are m x m and q x q diagonal matrices of the eigenvalue A, repeated rn and 
q times, respectively, i.e., 
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With the following simplifying definitions as in equations (31) 
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equation (Bl) becomes 
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where u(i) and y(i) are m x 1 and q x 1 input and output vectors, respectively. The above 
equation is in a linear form with the unknown parameters in the matrices a, J. and D with 
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Iii terms of the prescribed eigenvalues, has the following structure: 
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which has the same general structure as except the block matrices are of dimensions q x q. 
The mp x 1 input history vector u(i — p) and the qp x 1 output history vector y[i —p) are defined 
the same way as in equations (37) except that the input u(i) is of dimensions m x 1, and the 
output y(i) is of dimensions qx 1. Equation (B8) is in linear form; the parameter matrix 7 can be 
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solved for directly from input-output data. For on-line computation, the recursive least-squares 
solution to the parameter matrix 7 is given as 
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Finally, the actual system Markov parameters can then be recovered from the reconstructed 
observer Markov parameters according to equation ( 20 ) as 
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MIMO Complex Eigenvalue Assignment 

The complex eigenvalue assignment for the multiple-input, multiple-output case can be derived 
by setting A — T A C T , where A r is given as in equation (41 ). The prescribed complex conjugate 
pairs of eigenvalues are denoted A,- = ±ju t (t = 1,2,..., n/2). Using the same notation for 
vectors formed by the columns and rows of C* and /?*, respectively, the products in equation (Bl) 
may be expressed as 
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Similar definitions apply for a\^ and simply by replacing m by q. Equation (Bl) now 
becomes 
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Ill terms of the prescribed complex eigenvalues, %. m has the following structure: 



The structure for „ is similar. The recursive solution to equation (B18) can be obtained by 
replacing 7 by %, and T(i- 1) by rV(i- 1) in equations (B13). The observer Markov parameters 
and the actual Markov parameters can then be computed as 


Y r = CA T B = CA T B' -M 

= Qj&a §c ± 2 = YP Y { t 2) (B23) 
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Y T = Y? + Y, Y T } y T -i- 1 +Y [ t ] D 


v(2) ; 


i = 0 


— RrAr/m + 0c Ar^/^r-i -1 + Aty/£^ 


(H24) 


MIMO Mixed Real and Complex Eigenvalue Assignment 

Among n prescribed eigenvalues, let n r denote the number of prescribed real eigenvalues 
A i (i = 1, 2, .... n r ) and n c the number of prescribed complex eigenvalues a, ± JUtj 
(i = 1, 2, . . . , n f; /2). Then write A = T _1 A m T, and A m as in equation (Al), and define 



tim = 

& 

fir 

A (r) - 

Am.m — 

\M r A (r)/ 
Am Acju 

J 



j 




\ir) 



( r )' 


O' ' 


where 


'’(i) 6 (i) 




(«r) 


*7 

(«r) 


(n r +l) (n r + l) (»r“ 


2)^h 




— r* , 6 + / 

(h) <w- 



(B25) 
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a ( t) - 


A (r) A (r) ... A ,rj 

m.m — 2 ,rn — > 


A (r) - 
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\ {r) - 
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\l» dr) dr) 
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A (r) - 

Ar.q ~ 


T (r) . (r) fr) fr) 

-Id/ -1// — 2// — 2d/ 


n 7' 


_(r) 


Jr) 
/2d/ — n r /'2,q 


jm __ , (o) _ n 

’ —Lq ~ 7 — j.r/ “ 




Equation (Bl) may be expressed as 


P-1 P~1 

2/('*) = Rm ^ ^ — T — 1) + 0 m ^ ^ — 7” — 1) + Z?w(?) 

T=() 7—0 


7 mfm(i — 1 ) 


(B27) 
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Oth (J 1 ) 

7 m = QLm P m D r w(i ~ 1) = Vm{i ~ 1) 


(B28) 


<Pm (i ~ 1) = ^2 A mlnU{i - T - 1) - % nrn u(l - p) 


<£m(i - 1) = 5Z ~r-l) = %n^y(i ~ p) 

t= 0 J 

The matrix m includes elements formed from both real and complex prescribed eigenvalues: 


where 


(B30) 
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— n r /2,rn — n ( f‘2jn 


_( 1 ) j 

T ,, / 777 X 777 
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— 1,77! U 77!X77I 
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—2 777 •* X 777 


“L u mxm 


The diagonal matrices A- 7 ^, in m m are °f dimensions m x ra, i = 1, 2 , .... n r , 

or n r /2, and r = 1, 2, . . p — 1. Similar structures apply for which is composed 

of the matrices A - ^ u/j^ of dimensions q x q instead. The recursive solution to the 
parameter matrix 7 m is obvious. The observer Markov parameters and the actual system Markov 
parameters are simply 


\ r — \V ) Q XW/ 

r L*mAm.m Pm Am// 


F< 1) y<?> 


Y T — <* m A m . 7W + /3m { ^ + A 


(B31) 

(B32) 


MIMO Deadbeat Eigenvalue Assignment 

In the deadbeat case, all eigenvalues of the observer are placed at the origin. The 
corresponding Markov parameters will vanish identically after n time steps. In other words, 
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Y T = 0 for r = n, n + 1, n + 2, .... Let M ( i denote the deadbeat observer gain for the multiple- 
input multiple-output case. The input-output description is given in terms of the observer 
Markov parameters as 


n— 1 71 — i 

y(*) = y ( CATB d) u (* - r - !) - y ( CA T M d ) y(i - r - 1) + Du{i) 

r = 0 r=0 

n—l n— 1 

= y Y - T - 1) - Y r 2) y(i - T - 1) + Du(i) 

r=0 r— 0 

= - n) + fjdy ji - n) + Du{i) 

= 7uT«f(* - 1) (B33) 


where in the above equation B' = B + MjD, A = A + M ( {C . and 


7 d = 


Pd D 


r rf (* - 1) 


u(z — n) 
?/(* - n) 
u(i) 


(B34) 


a,/ = 


yW yW 
r 0 1 1 


Pd = 


y (2) v (^) 


y 1 



y( 2) 

1 n—l 




> 


(B35) 


The nm x 1 input history vector u(i — n) and the nq x 1 output history vector y(i — n) are 
defined as in equations (63) , except u(i) and y(i) are now m x 1 and g x 1 vectors, respectively. 
Note that in the deadbeat scheme, the observer Markov parameters are solved directly from 
input-output data, and the actual system Markov parameters are then recovered simply as in 
equation (20). 


41 



Appendix C 

The Mini-Mast Truss Structure 


A model obtained by finite element analysis of the Mini-Mast truss structure (ref. 21) is 
used as an example to illustrate the identification algorithms developed in this paper. The 
mathematical model has the first two bending inodes, with practically the same frequencies 
(0.8 Hz); the first torsional mode (4.3 Hz); and the second two bending modes (6.1 Hz), again 
with practically the same frequencies. The model considers two inputs and two outputs. The 
inputs are two torque wheels for the x arid y axes, and the outputs are two displacement sensors 
mounted at the top of the structure as shown in figure Cl. The system frequencies and the 
associated damping factors expressed as the real parts of the eigenvalues are listed in table Cl. 


Table Cl. Damping and Frequencies of 
Truss Structure 


Mode 

Damping factor 

Frequency, Hz 

i 

0.09 

0.80 

2 

.09 

.80 

3 

.33 

4.36 

4 

.38 

6.10 

5 

.39 

6.16 


The continuous-time system matrices are listed here. For ease of presentation, the matrices are 
subdivided and given below: 





where 



-8.918 x 10~ 2 

-1.330 x 10 -1 

-5.035 

4.756 x 10“ 5 

9.106 x 10 “ 4 

1.303 x 10~ 4 

-8.912 x 10~ 2 

-1.474 x 10~ 4 

5.032 

1.309 x 10" 2 

5.035 

1.540 x 10“ ‘ 

-9.212 x 10“ 2 

-1.293 x 10 4 

-1.403 x lO -3 

-4.100 x 10~ 5 

-5.032 

1.335 x 10~‘ 

-9.205 x 10^ 2 

-1.540 x lO -2 

-3.238 x 10~ 3 

-2.093 x 10~ 3 

3.540 x 10“ 3 

7.388 x 10~ 3 

-3.251 x 10 -1 

4.008 x 10~ 3 

-7.596 x HT 3 

-4.048 x 10 -3 

2.748 x 10" 3 

27.420 

2.468 x 10“ 2 

-9.535 x 10~ 2 

-2.691 x 10 ~ 2 

-1.040 x 10 _1 

1.546 x 10“ 3 

-9.585 x 10 -2 

-2.514 x 10~ 2 

1.043 x 10' 1 

-2.748 x 10~ 2 

3.791 x 10~ 3 

2.660 x 10~ 2 

-1.015 x 10 _1 

-2.617 x lO' 2 

-9.974 x 10- 2 

-3.283 x 10“ 3 

-1.020 x 10' 1 

-2.627 x 10“ 2 

1.005 x lO^ 1 

-2.567 x 10“ 2 

1.491 x 10 — 2 
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-1.549 x 10 -3 

5.498 x 10 3 

-1.999 x 10“ 2 

-9.892 x 10“ 3 

3.740 x 10~ 2 

1.527 x 10“ 2 

-2.110 x 10 2 

-6.062 x 10~ 3 

3.859 x 10~ 2 

1.079 x 10“ 2 

1.214 x 10~ 3 

1.010 x 10“ 2 

-3.766 x 10“ 2 

-5.495 x 10" 3 

2.055 x 10“ 2 

-1.412 x 1(T 2 

3.884 x 10-2 

1.075 x 10“ 2 

-2.179 x 10“ 2 

-6.517 x 10“ 3 

-27.420 

-9.447 x 10“ 3 

1.774 x 10“ 2 

1.114 x 10“ 2 

-2.519 x 10“ 2 

-3.330 x lO" 1 

-9.884 x 10“ 3 

2.238 x 10“ 2 

1.117 x 10~ 2 

-2.125 x 10 -2 

2.839 x 10 -3 

-3.763 x 10“' 

5.972 x 10 _1 

38.364 

-1.010 x 10" 1 

-1.320 x 10~ 2 

-5.956 x 10 _1 

-3.790 x 10“ 4 

-4.656 x 10“ 2 

38.660 

-2.834 x 10 -3 

-38.364 

4.638 x 10“ 2 

-3.912 x 10 1 

-5.969 x 10“ 1 

-2.657 x 10“ 3 

l.on x ur 1 

-38.660 

5.986 x 10 _1 

3.943 x lO' 1 


2.345 x 10 -3 

-1.996 x 10 '■ 

-2.101 x 10 -3 

-2.360 x 10 ; 

-2.349 x 10 -3 

1.999 x 10" 3 

-2.015 x 10 -3 

-2.364 x 10 '' 

-1.052 x 10 -4 

-2.488 x 10 _< 

1.107 x lO" 1 

2.455 x lO- 1 

1.667 x 10 3 

9.519 x 10 1 

-9.095 x 10“ 4 

1.554 x 10“ 3 

1.630 x 10“ 3 

9.180 x 10 4 

-8.917 x 10 1 

1.509 x 10 3 


and 


where 


= Ci 

C '2 

0.000 

0.000 

0.000 

0.000 



1.119 x 10" 2 

4.016 x 

10“ 3 

1.122 x 10 -2 

-4.025 x 10 3 

-9.167 x 10 

-9.114 x 10“ 3 

7.620 x 

10 -3 

-9.136 x 10~ 3 

-7.639 x 10 3 

-9.311 x 10 



-9.177 x 10 -3 

-4.321 x 10~ 4 

-2.448 x 10" 3 

4.669 x 10 4 

2.393 x 10 -3 

-9.326 x 10“ 3 

-2.427 x 10 -3 

1.965 x 10- 3 

2.423 x 10 3 

-1.990 x 10“ : 


In the numerical examples, the system model is discretized at a sampling frequency of *33.3 Hz. 
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